##############################
##Programme Laden
#############################
 
### Verzeichnist mit Modulen
    MODULE.DIR <- "C:/Documents and Settings/Stefan_Hoderlein/Desktop/"

### Einbinden der Funktionen
    source( paste(MODULE.DIR, "Gdata1.R", sep="") )
     source( paste(MODULE.DIR, "gen.af.R", sep="") )
 
 
 #########################################################################
 ## Daten aufbereiten                                                     
 #########################################################################
 filename <- "C:/Documents and Settings/Stefan_Hoderlein/Desktop/newGroups40.csv"
 
 A1 <- read.csv(filename , sep=";", dec=",")

##########HH Char
Zett <- cbind(A1$agehd1,A1$washmach,A1$ncars,A1$ifshhtyp)

Zschl <- gen.af(af=Zett,no=1)

##########Control Function

Endodat <- cbind(A1$totexp,A1$ln.hhincome,Zschl)
Endopos <- Endodat[,-1]
H1 <- c(15,0.5,0.3)
#Vau <-rep(NA,length(A1[,1]))

Endomod <- lm(Endodat[,1]~ Endodat[,2] + Endodat[,3])

summary(Endomod)

Vau <-Endomod$res

#Vau<-rnorm(mean=0,sd=1,n=length(A1[,1]))

summary(Vau)

#summary(A1)
Data <- cbind(A1$year,A1$FOOD/A1$relFOOD*100,A1$totexp,A1$relFOOD/100,Zschl,Vau)

Datap <- cbind(A1$year,A1$FOOD/A1$relFOOD*100,A1$HOUSING/A1$relHOUSING*100,A1$FUEL_LIG/A1$relFUEL_LIG*100,A1$relFOOD/100,A1$relHOUSING/100,A1$relFUEL_LIG/100)

#Data <- cbind(A1$year,A1$FOOD/A1$relFOOD*100,A1$totexp,A1$relFOOD/100)

colnames(Data)<-c("Year","Food","TotExp","Fprice","Princomp","controlF")
summary(Data)

Preis <-Data[,4]

for (j in 1:length(Data[,1]))

    {
    
    for (i in 1:27)

        {
        if (Data[j,1]==(1973+i))    
        Data[j,4] <- mean(Preis[Data[,1]==(1973+i)])
        }
    
    }

    
#################Group data so that we have 9 different scenarios

Daten1 <- matrix(Data[(Data[,1]==1974)|(Data[,1]==1975)|(Data[,1]==1976)],ncol=6)
Daten2 <- matrix(Data[(Data[,1]==1977)|(Data[,1]==1978)|(Data[,1]==1979)],ncol=6)
Daten3 <- matrix(Data[(Data[,1]==1980)|(Data[,1]==1981)|(Data[,1]==1982)],ncol=6)
Daten4 <- matrix(Data[(Data[,1]==1983)|(Data[,1]==1984)|(Data[,1]==1985)],ncol=6)
Daten5 <- matrix(Data[(Data[,1]==1986)|(Data[,1]==1987)|(Data[,1]==1988)],ncol=6)
Daten6 <- matrix(Data[(Data[,1]==1989)|(Data[,1]==1990)|(Data[,1]==1991)],ncol=6)
Daten7 <- matrix(Data[(Data[,1]==1992)|(Data[,1]==1993)|(Data[,1]==1994)],ncol=6)
Daten8 <- matrix(Data[(Data[,1]==1995)|(Data[,1]==1996)|(Data[,1]==1997)],ncol=6)
Daten9 <- matrix(Data[(Data[,1]==1998)|(Data[,1]==1999)|(Data[,1]==2000)],ncol=6)


Daten1[,4] <- rep(mean(Daten1[,4][Daten1[,1]==1975]),length(Daten1[,1]))
Daten2[,4] <- rep(mean(Daten2[,4][Daten2[,1]==1978]),length(Daten2[,1]))
Daten3[,4] <- rep(mean(Daten3[,4][Daten3[,1]==1981]),length(Daten3[,1]))
Daten4[,4] <- rep(mean(Daten4[,4][Daten4[,1]==1984]),length(Daten4[,1]))
Daten5[,4] <- rep(mean(Daten5[,4][Daten5[,1]==1987]),length(Daten5[,1]))
Daten6[,4] <- rep(mean(Daten6[,4][Daten6[,1]==1990]),length(Daten6[,1]))
Daten7[,4] <- rep(mean(Daten7[,4][Daten7[,1]==1993]),length(Daten7[,1]))
Daten8[,4] <- rep(mean(Daten8[,4][Daten8[,1]==1996]),length(Daten8[,1]))
Daten9[,4] <- rep(mean(Daten9[,4][Daten9[,1]==1999]),length(Daten9[,1]))

summary(Daten1)

##################calculate price change

deltap1 <-  Daten2[1,4] - Daten1[1,4]
deltap2 <-  Daten3[1,4] - Daten2[1,4]                                                                                                                                                                                                                                                                                                                          
deltap3 <-  Daten4[1,4] - Daten3[1,4]
deltap4 <-  Daten5[1,4] - Daten4[1,4]
deltap5 <-  Daten6[1,4] - Daten5[1,4]
deltap6 <-  Daten7[1,4] - Daten6[1,4]                                                                                                                                                                                                                                                                                                                          
deltap7 <-  Daten8[1,4] - Daten7[1,4]
deltap8 <-  Daten9[1,4] - Daten8[1,4] 

print(c(deltap1,deltap2,deltap3,deltap4,deltap5,deltap6,deltap7,deltap8))

#########prepare the variables

Q2maldeltap2 <- Daten2[,2]*deltap2
Q3maldeltap2 <- Daten3[,2]*deltap2

###############This is the Left hand side

print(summary(Q2maldeltap2))

plot(density(Q2maldeltap2),xlab="Density of deltap*Q")

B1<- sort(Daten3[,3],index.return=T)
B2<-B1$ix

IndeX <- rep(NA,10)

for (s in 0:9)
    
    {IndeX[s+1] <- round(length(Daten3[,3])*(5+10*s)/100)}
    
Index <- B2[IndeX]
Quant3 <- Daten3[Index,3]
Quant3


################# End Quant t+1

B3<- sort(Daten2[,3],index.return=T)
B4<-B3$ix

IndeX <- rep(NA,10)

for (s in 0:9)
    
    {IndeX[s+1] <- round(length(Daten2[,3])*(5+10*s)/100)}
    
Index2 <- B4[IndeX]
Quant2 <- Daten2[Index2,3]
Quant2

############### End Quant t

Quant2 <- 0.5*(Quant3 + Quant2)



############These are points on Y2 vector


#Positions <- cbind(rep(Quant2,10),c(rep(Quant3[1],10),rep(Quant3[2],10),rep(Quant3[3],10),rep(Quant3[4],10),rep(Quant3[5],10),rep(Quant3[6],10),rep(Quant3[7],10),rep(Quant3[8],10),rep(Quant3[9],10),rep(Quant3[10],10)))

Quant4 <- Quant2 - 1.1



for (s in 1:9)

    {Quant4<-c(Quant4,Quant2 - 1.1 + 0.1*s)}
  
    Positions <- cbind(rep(Quant2,10),Quant4)
 

Diff <- Positions[,2] - Positions[,1]
Position <- cbind(Positions,Diff,rep(0,100),rep(-5,100))
Position
#H1<-c(20,1.5,20)


#Y2minY1<-matrix(NA,ncol=10,nrow=10) 
#for (i in 1:10)
#{Y2 <- Quant3[i]
#    for (j in 1:10)
#       {Y2minY1[j,i] <- Y2 - Quant2[j]
#        Positions[,3] <- Y2 - Quant2[j]}
#    }

#print(summary(Y2minY1)
print(summary(Q2maldeltap2))

##############These are Differences between Y2 and Y1

#########Now the probs

P2minP1 <- matrix(NA,nrow=100,ncol=4) 
BootSE  <- matrix(NA,nrow=100,ncol=4)
#P2minP1<-matrix(NA,nrow=100,ncol=6)
refine<-1


for (k in 1:100)

   {    print(k)
        Ind1 <- sign(Position[k,3]-Q2maldeltap2)
        Ind1[Ind1<0] <- 0  
        Dat1 <- cbind(Ind1,Daten2[,3])
        
        Erg1 <- Gdata(data=Dat1,h=H1,x0=c(Position[k,1],0,0),noend=1) 
        Erg1 <- Erg1[1]
            
        Ind2 <- sign(Position[k,3]-Q3maldeltap2)
        Ind2[Ind2<0] <- 0 
        Dat2 <- cbind(Ind2,Daten3[,3])
        
        Erg2 <- Gdata(data=Dat2,h=H1,x0=c(Position[k,2],0,0),noend=1) 
        Erg2 <- Erg2[1]
       
        
         if (refine==1)
        
        {
       P2minP1[k,] <- c(Erg1,Erg2,max(Erg1-Erg2,0),min(Erg1,1-Erg2))
        }
        
        else
        {
          P2minP1[k,] <- c(Erg1,Erg2,max(Erg1-Erg2,0),Erg1*(1-Erg2))
        }
        ######################Bootstrap 
        
        Bootse <- matrix(NA,nrow=101,ncol=3)
           
           for (l in 1:101) 
           {
            #print(l)
            x<-c(1:length(Dat1[,1]))
            IndX1<-sample(x,replace=T)
            Dat1boot<-Dat1[IndX1,]
            ERG1 <- Gdata(data=Dat1boot,h=H1,x0=c(Position[k,1],0,0),noend=1) 
            ERG1 <- ERG1[1]
                
            y<-c(1:length(Dat2[,1]))
            IndX2<-sample(y,replace=T)
            Dat2boot<-Dat2[IndX2,]
            ERG2 <- Gdata(data=Dat2boot,h=H1,x0=c(Position[k,1],0,0),noend=1) 
            ERG2 <- ERG2[1]  
            #############lower bound
            
            
            ERGrest <- 0.5*(ERG1+ERG2)
            Ystern <- c(rep(1,round(ERGrest*100)),rep(0,round((1-ERGrest)*100)))
            Ystern1 <- sample (Ystern,size=length(Ystern), replace=F)
            Dat3boot <- cbind(Ystern1,Dat1boot[,-1])
            Dat4boot <- cbind(Ystern1,Dat2boot[,-1])
            ERG3 <- Gdata(data=Dat3boot,h=H1,x0=c(Position[k,1],0,0),noend=1) 
            ERG3 <- ERG3[1]
                
            ERG4 <- Gdata(data=Dat4boot,h=H1,x0=c(Position[k,1],0,0),noend=1) 
            ERG4 <- ERG4[1] 
            
            Bootse[l,] <- c(ERG1-ERG2-(Erg1-Erg2),ERG1-Erg1,Erg2-ERG2)                
            phin <-(length(Dat1[,1])*H1[1]*H1[2]*H1[3])^0.5
            Bootse[l,] <- phin*Bootse[l,]
            }
            
            u5  <-order(Bootse[,1])
            BU  <- Bootse[u5,1]
            BUU <- Erg1-Erg2 - (1/phin)*BU[95]
            BUO <- Erg1-Erg2 - (1/phin)*BU[3]
        
            O5  <-order(Bootse[,2])
            BO  <- Bootse[O5,2]
            TestA1<-Erg1-(1/phin)*BO[90]
            TestA2<-Erg1-(1/phin)*BO[10]
            #BOU <- Erg1-(1/phin)*BO[98]
            #BOO <- Erg1-(1/phin)*BO[3]
            
            O52  <-order(Bootse[,3])
            BO2  <- Bootse[O52,3]
            TestB1<- 1-Erg2-(1/phin)*BO2[90]
            TestB2<- 1-Erg2-(1/phin)*BO2[10]
            #BOU2 <- BO2[2]
            #BOO2 <- BO2[98]
         
        if (refine==1)
        
        {
         if (TestA1>TestB2)
            (
            BootSE[k,] <- c(max(BUU,0),BUO,max(1-Erg2-(1/phin)*BO2[97],0),1-Erg2-(1/phin)*BO2[3])
            )
            else  if (TestB1>TestA2)
                {    
                BootSE[k,] <- c(max(BUU,0),BUO,max(Erg1-(1/phin)*BO[97],0),Erg1-(1/phin)*BO[3])
                }
                else
                    {
                    BootSE[k,] <- c(max(BUU,0),BUO,max(min(Erg1,1-Erg2)-(1/phin)*BO[98],0),min(Erg1,1-Erg2)-(1/phin)*BO[2])
                    }
         }
         
         else {BootSE[k,] <- c(max(BUU,0),BUO,max(Erg1*(1-Erg2)-(1/phin)*BO[98],0),Erg1*(1-Erg2)-(1/phin)*BO[2])}
   }
    
ERGEBNIS<-cbind(Position,P2minP1,BootSE)
ERGEBNIS
matrix(ERGEBNIS[ERGEBNIS[,8]>0.01|ERGEBNIS[,9]>0.01],ncol=13)

plot(density(ERGEBNIS[,9][ERGEBNIS[,9]>0.01]),type="l",ylab="Rel.Frequency Upper Bound",xlab="Point Estimate of Upper Bound at Varying Values of Regressors", main="Distribution of Upper Bound and Associated Lower Conf. Int.")
lines(density(ERGEBNIS[,12][ERGEBNIS[,9]>0.01]),lty=2)

